
% Here we want to find R (June 11th 2021)

function R = findR(Q, I, lambda, Nmax)
R =zeros(I, Nmax+1);

% First, calculate the (1,1) value.

% R(1,1) = Q(1,1);
for m=0 : Nmax
    R(1,1) = R(1,1)+(1-poisscdf (m+1, lambda) ) * Q(1,m+1);
    % NOTICE! Here m should count 1 out.
end

% Then for the rest of row 1
for m=1 : Nmax
	for k =1:Nmax       
        R(1,k+1)=R(1,k+1)+poisspdf (m-k,lambda) * Q(1,m+1);
	end
end

% For row 2 to n
for i=2:I
    % For first column
	for m=0:Nmax     % Here k2 infers ktilt
        for k2=0:Nmax
            R(i,1)=R(i,1)+(1-poisscdf (m+ k2+1, lambda)) * Q(i,m+1) * R(i-1,k2+1); 
        end
    end
    
  % Fir column 2 to Nmax+1
	for k=1: Nmax
        for m=0 : Nmax
            for k2=0:Nmax       % Here k2 infers ktilt
                R(i,k+1)=R(i,k+1) + poisspdf(m+k2-k,lambda)*Q(i,m+1)*R(i-1,k2+1);
            end
        end
	end
end
